# import packages
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import sys
from pandas._libs.lib import is_integer

## Linear Estimator
def regEstimate(dataset, pbvarb = 'predicted_prob_black', outcome = 'audited', wvarb = None):
	if wvarb is not None: 
		model =  smf.wls(outcome + ' ~ ' + pbvarb, dataset, weights=dataset[wvarb]).fit(cov_type='HC1')
		coef= model.params[pbvarb]
		se = model.bse[pbvarb]
		return coef, se, model
	else:
		model =  smf.ols(outcome + ' ~ ' + pbvarb, dataset).fit(cov_type='HC1')
		coef= model.params[pbvarb]
		se = model.bse[pbvarb]
		return coef, se, model

## Probabilistic Estimator 
def chenEstimate(dataset, pbvarb='pBlack', outcome='audited', wvarb=None):
	if wvarb is not None:
		return (dataset[pbvarb]*dataset[outcome]*dataset[wvarb]).sum()/(dataset[pbvarb]*dataset[wvarb]).sum()-((1-dataset[pbvarb])*dataset[outcome]*dataset[wvarb]).sum()/((1-dataset[pbvarb])*dataset[wvarb]).sum()
	else:
		return (dataset[pbvarb]*dataset[outcome]).sum()/(dataset[pbvarb]).sum()-((1-dataset[pbvarb])*dataset[outcome]).sum()/(1-dataset[pbvarb]).sum()

## functions for calculataxpayer_idg standard errors
def getWVar(values, weights):
	average = np.average(values, weights=weights)
	variance = np.average((values-average)**2, weights=weights)
	return variance

def getSEMultiplier(dataset, pbvarb='pBlack', wvarb=None):
	if wvarb is not None:
		return np.sqrt(getWVar(dataset[pbvarb], dataset[wvarb])/((dataset[pbvarb]*dataset[wvarb]).mean()*((1-dataset[pbvarb])*dataset[wvarb]).mean()))
	else:
		return np.sqrt(dataset[pbvarb].var()/(dataset[pbvarb].mean()*(1-dataset[pbvarb].mean())))

def getSEs(dataset, pbvarb='pBlack', outcome='audited', wvarb=None):
	seMultiplier=getSEMultiplier(dataset,pbvarb,wvarb=wvarb)
	seReg = regEstimate(dataset,pbvarb,outcome, wvarb=wvarb)[1]
	seChen = seReg*seMultiplier
	return seChen, seReg

def w_avg(df, values, weights):
	d = df[values]
	w = df[weights]
	return (d * w).sum()/w.sum()

# function for calculataxpayer_idg covariance conditions
def get_cond_cov(df, covarb, conditionvarb, outcomevar ,wgtvar):
        #print('startaxpayer_idg len of dataset: ' + str(df.shape[0]))
        df = df.dropna(subset = [covarb, conditionvarb, outcomevar])
        #print('len of data after dropping nulls: ' + str(df.shape[0]))
        conditions = df[conditionvarb].unique()
        totalwt=df[wgtvar].sum()
        #print('total weight is: ' + str(totalwt))
        estimates = []
        #print('conditioning varbs: ' + str(conditions))
        for cond in conditions:
                 estimate = {}
                 estimate['cond_val'] = cond
                 dat = df[df[conditionvarb] == cond]
                 estimate['p_cond'] = (dat[wgtvar]/totalwt).sum()
                 reg = smf.wls(outcomevar+' ~ '+covarb, data= dat, weights= dat[wgtvar]).fit(cov_type='HC1')
                 estimate['var'] = dat[covarb].var()
                 estimate['coef'] = reg.params[covarb]
                 estimates.append(estimate)
        estimates = pd.DataFrame(estimates)
        exp_cond_cov = (estimates['p_cond']*estimates['coef']*estimates['var']).sum()
        return exp_cond_cov

def getWVar(values, weights):
        average = np.average(values, weights=weights)
        variance = np.average((values-average)**2, weights=weights)
        return variance


## files to read: 

# individual2014 
# ncdata 
# ncweights 

keep_cols = ['taxpayer_id', 'predicted_prob_black', 'isEIC', 'filing_jointly', 'isM', 'pos_deps', 'aud_no_research_audits', 'adj_gross_inc']
individual2014 = pd.read_csv('/REDACTED/data/final/individualBISG2014_full_final.csv', usecols=keep_cols)
ncdata = pd.read_csv("/REDACTED/NC_Analysis_Dataset_2014_tplevel.csv")
ncweights = pd.read_csv("/REDACTED/HertzGraff/nc_reweights_v4_April.csv")

# clean datasets and merge together
ncdata = ncdata[ncdata.black_ind.notna()]
ncdata = ncdata[['taxpayer_id', 'black_ind']]
ncweights = ncweights[['taxpayer_id', 'uswgt']]

individual2014 = pd.merge(individual2014, ncdata, on = 'taxpayer_id', how = 'left')
len(individual2014[individual2014.black_ind.notnull()])
len(individual2014)
individual2014 = pd.merge(individual2014, ncweights, on = 'taxpayer_id', how = 'left')
len(individual2014[individual2014.uswgt.notnull()])
len(individual2014)

# clean merged dataset
individual2014['unitwt'] = 1
individual2014['p_black_rd'] = individual2014['predicted_prob_black'].round(2)
individual2014['p_black_rd'] = pd.to_numeric(individual2014['p_black_rd'])
individual2014['predicted_prob_black'] = pd.to_numeric(individual2014['predicted_prob_black'])
individual2014['black_ind'] = pd.to_numeric(individual2014['black_ind'])

######################
### EITC SUBGROUPS ###
######################

# define subgroups for recalibration
dfname = ['Single EITC', 'Single EITC wgt',
           'Joint EITC', 'Joint EITC wgt',
           'Single Male EITC', 'Single Male EITC wgt',
           'Single Female EITC', 'Single Female EITC wgt',
           'Single Male EITC with Deps', 'Single Male EITC with Deps wgt',
           'Single Male EITC without Deps', 'Single Male EITC without Deps wgt'] 
dflist = [(individual2014.isEIC == 1) & ~(individual2014.filing_jointly == 1),
           (individual2014.isEIC == 1) & ~(individual2014.filing_jointly == 1),
           (individual2014.isEIC == 1) & (individual2014.filing_jointly == 1),
           (individual2014.isEIC == 1) & (individual2014.filing_jointly == 1),
           (individual2014.isEIC == 1) & ~(individual2014.filing_jointly == 1) & (individual2014.isM == 1),
           (individual2014.isEIC == 1) & ~(individual2014.filing_jointly == 1) & (individual2014.isM == 1),
           (individual2014.isEIC == 1) & ~(individual2014.filing_jointly == 1) & ~(individual2014.isM == 1),
           (individual2014.isEIC == 1) & ~(individual2014.filing_jointly == 1) & ~(individual2014.isM == 1),
           (individual2014.isEIC == 1) & ~(individual2014.filing_jointly == 1) & (individual2014.isM == 1) & (individual2014.pos_deps == 1),
           (individual2014.isEIC == 1) & ~(individual2014.filing_jointly == 1) & (individual2014.isM == 1) & (individual2014.pos_deps == 1),
           (individual2014.isEIC == 1) & ~(individual2014.filing_jointly == 1) & (individual2014.isM == 1)  & ~(individual2014.pos_deps == 1),
           (individual2014.isEIC == 1) & ~(individual2014.filing_jointly == 1) & (individual2014.isM == 1)  & ~(individual2014.pos_deps == 1)]


# write out results
sys.stdout = open("/REDACTED/final_NC_cov_est_full_subgroups_align_V3.txt", "w")

# loop over subgroups and calculate recalibration results
i = 0
for dfs in dflist:
        print(dfname[i])
        i = i+1
        df = individual2014[dfs].copy()
        nc_df = df[(df.black_ind.notnull()) & (df.uswgt.notnull())]
        ## recalibrate prob_black
        B_b_model = sm.WLS.from_formula("black_ind ~ predicted_prob_black", data = nc_df, weights = nc_df.unitwt).fit(cov_type = "HC1")
        B_b_model_weight = sm.WLS.from_formula("black_ind ~ predicted_prob_black", weights = nc_df.uswgt, data = nc_df).fit(cov_type = "HC1")
        #### rho
        rho = B_b_model.params['predicted_prob_black']
        rho_weight = B_b_model_weight.params['predicted_prob_black']
        df['bhat'] = B_b_model.predict(df['predicted_prob_black'])
        df['epsilonhat'] = B_b_model.resid
        df['bhat_weight'] = B_b_model_weight.predict(df['predicted_prob_black'])
        df['epsilonhat_weight'] = B_b_model_weight.resid
        lin_disp = regEstimate(df, pbvarb = 'bhat', outcome = 'aud_no_research_audits', wvarb='unitwt')
        lin_disp_weight = regEstimate(df, pbvarb = 'bhat_weight', outcome = 'aud_no_research_audits', wvarb='uswgt')   
        df['etahat'] = lin_disp[2].resid
        df['etahat_weight'] = lin_disp_weight[2].resid
        print("done")
        #### cov(b, epsilon)
        df['b_eps_product'] =df['epsilonhat']*df['predicted_prob_black']
        df['b_eps_product_weight'] =df['epsilonhat_weight']*df['predicted_prob_black']
        cov_b_eps = w_avg(df, 'b_eps_product', 'unitwt') - w_avg(df, 'epsilonhat', 'unitwt')*w_avg(df, 'predicted_prob_black', 'unitwt')
        cov_b_eps_weight = w_avg(df, 'b_eps_product_weight', 'uswgt') -  w_avg(df, 'epsilonhat_weight', 'uswgt')* w_avg(df, 'predicted_prob_black', 'uswgt') 
        #### Covariance term Cov(E(eta|bhat), E(epsilon|bhat))
        df["bhat_bucket"] = df.bhat.round(2)
        mu_epsilon_dict = df.groupby('bhat_bucket').epsilonhat.mean().to_dict()
        mu_eta_dict = df.groupby('bhat_bucket').etahat.mean().to_dict()
        df['mu_epsilon'] = df.bhat_bucket.map(mu_epsilon_dict).astype(float)
        df['mu_eta'] = df.bhat_bucket.map(mu_eta_dict).astype(float)
        df['eta_epsilon_product'] = df['mu_epsilon'] * df['mu_eta']	
        overall_cov_term = w_avg(df, 'eta_epsilon_product', 'unitwt') - w_avg(df, 'mu_eta', 'unitwt')*w_avg(df, 'mu_epsilon', 'unitwt')
        df["bhat_bucket_weight"] = df.bhat_weight.round(2)
        mu_epsilon_dict_weight = df.groupby('bhat_bucket_weight').apply(w_avg, 'epsilonhat_weight', 'uswgt').to_dict()
        mu_eta_dict_weight = df.groupby('bhat_bucket_weight').apply(w_avg, 'etahat_weight', 'uswgt').to_dict()
        df['mu_epsilon_weight'] = df.bhat_bucket_weight.map(mu_epsilon_dict_weight).astype(float)
        df['mu_eta_weight'] = df.bhat_bucket_weight.map(mu_eta_dict_weight).astype(float)
        df['eta_epsilon_product_weight'] = df['mu_epsilon_weight'] * df['mu_eta_weight']	
        overall_cov_term_weight =  w_avg(df, 'eta_epsilon_product_weight', 'uswgt') - w_avg(df, 'mu_eta_weight', 'uswgt')*w_avg(df, 'mu_epsilon_weight', 'uswgt')
        ### E[cov(Y,b|B)]
        cov_Y_b_B = get_cond_cov(df, 'predicted_prob_black', 'black_ind', 'aud_no_research_audits' ,'unitwt')
        cov_Y_b_B_weight = get_cond_cov(df, 'predicted_prob_black', 'black_ind', 'aud_no_research_audits' ,'uswgt')
        ### E[cov(Y,b*|B)]
        cov_Y_bhat_B = get_cond_cov(df, 'bhat', 'black_ind', 'aud_no_research_audits' ,'unitwt')
        cov_Y_bhat_B_weight = get_cond_cov(df, 'bhat_weight', 'black_ind', 'aud_no_research_audits' ,'uswgt')
        ### E[cov(Y,B|b)]
        cov_Y_B_b = get_cond_cov(df, 'black_ind', 'p_black_rd', 'aud_no_research_audits' ,'unitwt')
        cov_Y_B_b_weight = get_cond_cov(df, 'black_ind', 'p_black_rd', 'aud_no_research_audits' ,'uswgt')
        ### E[cov(Y,B|b*)]
        cov_Y_B_bhat = get_cond_cov(df, 'black_ind', 'bhat_bucket', 'aud_no_research_audits' ,'unitwt')
        cov_Y_B_bhat_weight = get_cond_cov(df, 'black_ind', 'bhat_bucket_weight', 'aud_no_research_audits' ,'uswgt')
        ### Truth
        D = w_avg(df[df['black_ind'] == 1], 'aud_no_research_audits', 'unitwt') - w_avg(df[df['black_ind'] == 0], 'aud_no_research_audits', 'unitwt')
        D_weight = w_avg(df[df['black_ind'] == 1], 'aud_no_research_audits', 'uswgt') - w_avg(df[df['black_ind'] == 0], 'aud_no_research_audits', 'uswgt')
        ### D*_l
        D_l = lin_disp[2].params['bhat']
        D_l_weight = lin_disp_weight[2].params['bhat_weight']
        ### D*_p
        D_p = chenEstimate(df, pbvarb = 'bhat', outcome = 'aud_no_research_audits', wvarb= 'unitwt')
        D_p_weight = chenEstimate(df, pbvarb = 'bhat_weight', outcome = 'aud_no_research_audits', wvarb= 'uswgt')
        var_bhat = getWVar(df.bhat, df.unitwt)
        var_bhat_weight = getWVar(df.bhat, df.uswgt)
        var_B = getWVar(df.black_ind, df.unitwt)
        var_B_weight = getWVar(df.black_ind, df.uswgt)
        print("Unweighted")
        print("rho", rho)
        print("var_bhat", var_bhat)
        print("var_B", var_B)
        print("cov_b_eps", cov_b_eps)
        print("cov_Y_b_B", cov_Y_b_B)
        print("cov_Y_B_b", cov_Y_B_b)
        print("cov_Y_bhat_B", cov_Y_bhat_B)
        print("cov_Y_B_bhat", cov_Y_B_bhat)
        print("overall_cov_term", overall_cov_term)
        print("D_l", D_l)
        print("D", D)
        print("D_p", D_p)
        print("Weighted")
        print("rho", rho_weight)
        print("var_bhat", var_bhat_weight)
        print("var_B", var_B_weight)
        print("cov_b_eps", cov_b_eps_weight)
        print("cov_Y_b_B", cov_Y_b_B_weight)
        print("cov_Y_B_b", cov_Y_B_b_weight)
        print("cov_Y_bhat_B", cov_Y_bhat_B_weight)
        print("cov_Y_B_bhat", cov_Y_B_bhat_weight)
        print("overall_cov_term", overall_cov_term_weight)
        print("D_l", D_l_weight)
        print("D", D_weight)
        print("D_p", D_p_weight)
        print("\n\n\n")

sys.stdout = sys.__stdout__

###################
### INCOME BINS ###
###################

# subset to taxpayers with non-negative AGI and split into 20 AGI bins
dataBISG_noneg_agi = individual2014[individual2014['adj_gross_inc']>=0]
n_bin = 20
dataBISG_noneg_agi["bin"] = pd.qcut(dataBISG_noneg_agi["adj_gross_inc"], q = n_bin, labels = list(range(1, n_bin+1)))
dataBISG_noneg_agi.bin.value_counts()
dataBISG_noneg_agi["bin"] = dataBISG_noneg_agi.bin.astype(float).astype('Int64')
bin_20_mean = dataBISG_noneg_agi.groupby('bin')['adj_gross_inc'].mean()

# define income bins for recalibration
dfname = ["Bin 1",
          "Bin 2", 
          "Bin 3",
          "Bin 4",
          "Bin 5",
          "Bin 6",
          "Bin 7",
          "Bin 8",
          "Bin 9",
          "Bin 10",
          "Bin 11",
          "Bin 12",
          "Bin 13",
          "Bin 14",
          "Bin 15",
          "Bin 16",
          "Bin 17",
          "Bin 18",
          "Bin 19",
          "Bin 20"] 

dflist = [dataBISG_noneg_agi.bin == 1,
          dataBISG_noneg_agi.bin == 2,
          dataBISG_noneg_agi.bin == 3,
          dataBISG_noneg_agi.bin == 4,
          dataBISG_noneg_agi.bin == 5,
          dataBISG_noneg_agi.bin == 6,
          dataBISG_noneg_agi.bin == 7,
          dataBISG_noneg_agi.bin == 8,
          dataBISG_noneg_agi.bin == 9,
          dataBISG_noneg_agi.bin == 10,
          dataBISG_noneg_agi.bin == 11,
          dataBISG_noneg_agi.bin == 12,
          dataBISG_noneg_agi.bin == 13,
          dataBISG_noneg_agi.bin == 14,
          dataBISG_noneg_agi.bin == 15,
          dataBISG_noneg_agi.bin == 16,
          dataBISG_noneg_agi.bin == 17,
          dataBISG_noneg_agi.bin == 18,
          dataBISG_noneg_agi.bin == 19,
          dataBISG_noneg_agi.bin == 20]

# write out results
sys.stdout = open("/REDACTED/final_NC_cov_est_full_bins_align_V3.txt", "w")

# loop over income bins and calculate recalibration results
i = 0
for dfs in dflist:
        print(dfname[i])
        i = i+1
        df = dataBISG_noneg_agi[dfs].copy()
        nc_df = df[(df.black_ind.notnull()) & (df.uswgt.notnull())]
        ## recalibrate prob_black
        B_b_model = sm.WLS.from_formula("black_ind ~ predicted_prob_black", data = nc_df, weights = nc_df.unitwt).fit(cov_type = "HC1")
        B_b_model_weight = sm.WLS.from_formula("black_ind ~ predicted_prob_black", weights = nc_df.uswgt, data = nc_df).fit(cov_type = "HC1")
        #### rho
        rho = B_b_model.params['predicted_prob_black']
        rho_weight = B_b_model_weight.params['predicted_prob_black']
        df['bhat'] = B_b_model.predict(df['predicted_prob_black'])
        df['epsilonhat'] = B_b_model.resid
        df['bhat_weight'] = B_b_model_weight.predict(df['predicted_prob_black'])
        df['epsilonhat_weight'] = B_b_model_weight.resid
        lin_disp = regEstimate(df, pbvarb = 'bhat', outcome = 'aud_no_research_audits', wvarb='unitwt')
        lin_disp_weight = regEstimate(df, pbvarb = 'bhat_weight', outcome = 'aud_no_research_audits', wvarb='uswgt')   
        df['etahat'] = lin_disp[2].resid
        df['etahat_weight'] = lin_disp_weight[2].resid
        print("done")
        #### cov(b, epsilon)
        df['b_eps_product'] =df['epsilonhat']*df['predicted_prob_black']
        df['b_eps_product_weight'] =df['epsilonhat_weight']*df['predicted_prob_black']
        cov_b_eps = w_avg(df, 'b_eps_product', 'unitwt') - w_avg(df, 'epsilonhat', 'unitwt')*w_avg(df, 'predicted_prob_black', 'unitwt')
        cov_b_eps_weight = w_avg(df, 'b_eps_product_weight', 'uswgt') -  w_avg(df, 'epsilonhat_weight', 'uswgt')* w_avg(df, 'predicted_prob_black', 'uswgt') 
        #### Covariance term Cov(E(eta|bhat), E(epsilon|bhat))
        df["bhat_bucket"] = df.bhat.round(2)
        mu_epsilon_dict = df.groupby('bhat_bucket').epsilonhat.mean().to_dict()
        mu_eta_dict = df.groupby('bhat_bucket').etahat.mean().to_dict()
        df['mu_epsilon'] = df.bhat_bucket.map(mu_epsilon_dict).astype(float)
        df['mu_eta'] = df.bhat_bucket.map(mu_eta_dict).astype(float)
        df['eta_epsilon_product'] = df['mu_epsilon'] * df['mu_eta']	
        overall_cov_term = w_avg(df, 'eta_epsilon_product', 'unitwt') - w_avg(df, 'mu_eta', 'unitwt')*w_avg(df, 'mu_epsilon', 'unitwt')
        df["bhat_bucket_weight"] = df.bhat_weight.round(2)
        mu_epsilon_dict_weight = df.groupby('bhat_bucket_weight').apply(w_avg, 'epsilonhat_weight', 'uswgt').to_dict()
        mu_eta_dict_weight = df.groupby('bhat_bucket_weight').apply(w_avg, 'etahat_weight', 'uswgt').to_dict()
        df['mu_epsilon_weight'] = df.bhat_bucket_weight.map(mu_epsilon_dict_weight).astype(float)
        df['mu_eta_weight'] = df.bhat_bucket_weight.map(mu_eta_dict_weight).astype(float)
        df['eta_epsilon_product_weight'] = df['mu_epsilon_weight'] * df['mu_eta_weight']	
        overall_cov_term_weight =  w_avg(df, 'eta_epsilon_product_weight', 'uswgt') - w_avg(df, 'mu_eta_weight', 'uswgt')*w_avg(df, 'mu_epsilon_weight', 'uswgt')
        ### E[cov(Y,b|B)]
        cov_Y_b_B = get_cond_cov(df, 'predicted_prob_black', 'black_ind', 'aud_no_research_audits' ,'unitwt')
        cov_Y_b_B_weight = get_cond_cov(df, 'predicted_prob_black', 'black_ind', 'aud_no_research_audits' ,'uswgt')
        ### E[cov(Y,b*|B)]
        cov_Y_bhat_B = get_cond_cov(df, 'bhat', 'black_ind', 'aud_no_research_audits' ,'unitwt')
        cov_Y_bhat_B_weight = get_cond_cov(df, 'bhat_weight', 'black_ind', 'aud_no_research_audits' ,'uswgt')
        ### E[cov(Y,B|b)]
        cov_Y_B_b = get_cond_cov(df, 'black_ind', 'p_black_rd', 'aud_no_research_audits' ,'unitwt')
        cov_Y_B_b_weight = get_cond_cov(df, 'black_ind', 'p_black_rd', 'aud_no_research_audits' ,'uswgt')
        ### E[cov(Y,B|b*)]
        cov_Y_B_bhat = get_cond_cov(df, 'black_ind', 'bhat_bucket', 'aud_no_research_audits' ,'unitwt')
        cov_Y_B_bhat_weight = get_cond_cov(df, 'black_ind', 'bhat_bucket_weight', 'aud_no_research_audits' ,'uswgt')
        ### Truth
        D = w_avg(df[df['black_ind'] == 1], 'aud_no_research_audits', 'unitwt') - w_avg(df[df['black_ind'] == 0], 'aud_no_research_audits', 'unitwt')
        D_weight = w_avg(df[df['black_ind'] == 1], 'aud_no_research_audits', 'uswgt') - w_avg(df[df['black_ind'] == 0], 'aud_no_research_audits', 'uswgt')
        ### D*_l
        D_l = lin_disp[2].params['bhat']
        D_l_weight = lin_disp_weight[2].params['bhat_weight']
        ### D*_p
        D_p = chenEstimate(df, pbvarb = 'bhat', outcome = 'aud_no_research_audits', wvarb= 'unitwt')
        D_p_weight = chenEstimate(df, pbvarb = 'bhat_weight', outcome = 'aud_no_research_audits', wvarb= 'uswgt')
        var_bhat = getWVar(df.bhat, df.unitwt)
        var_bhat_weight = getWVar(df.bhat, df.uswgt)
        var_B = getWVar(df.black_ind, df.unitwt)
        var_B_weight = getWVar(df.black_ind, df.uswgt)
        print("Unweighted")
        print("rho", rho)
        print("var_bhat", var_bhat)
        print("var_B", var_B)
        print("cov_b_eps", cov_b_eps)
        print("cov_Y_b_B", cov_Y_b_B)
        print("cov_Y_B_b", cov_Y_B_b)
        print("cov_Y_bhat_B", cov_Y_bhat_B)
        print("cov_Y_B_bhat", cov_Y_B_bhat)
        print("overall_cov_term", overall_cov_term)
        print("D_l", D_l)
        print("D", D)
        print("D_p", D_p)
        print("bin_20_mean", bin_20_mean[i])
        print("Weighted")
        print("rho", rho_weight)
        print("var_bhat", var_bhat_weight)
        print("var_B", var_B_weight)
        print("cov_b_eps", cov_b_eps_weight)
        print("cov_Y_b_B", cov_Y_b_B_weight)
        print("cov_Y_B_b", cov_Y_B_b_weight)
        print("cov_Y_bhat_B", cov_Y_bhat_B_weight)
        print("cov_Y_B_bhat", cov_Y_B_bhat_weight)
        print("overall_cov_term", overall_cov_term_weight)
        print("D_l", D_l_weight)
        print("D", D_weight)
        print("D_p", D_p_weight)
        print("bin_20_mean", bin_20_mean[i])
        print("\n\n\n")


sys.stdout = sys.__stdout__